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EXECUTIVE  SUMMARY 


The  goal  of  this  project  is  to  develop  a  unified  computational  framework  for  multi-scale 
simulations  of  weakly-ionized  non-equilibrium  plasmas.  This  Final  Report  for  Phase  I  covers 
the  performance  period  from  September  1,  2007  through  June  30,  2008.  The  main  results  of 
the  work  can  be  summarized  as  follows 

•  We  have  investigated  methods  of  Renormalization  Group  and  Slow  Invariant 
Manifolds  in  application  to  plasma  modeling.  Using  the  point  symmetry  group  for  the 
Vlasov-Maxwell  system,  we  obtained  exact  solutions  for  the  dynamics  of  collisionless 
plasma.  The  interpretation  of  the  solutions  in  terms  of  invariants  of  the  group 
transformation  and  invariant  manifolds  was  achieved  -  the  particle  distribution 
functions  were  related  to  invariants  of  the  group  transform.  The  road  towards  obtaining 
approximate  symmetries  due  to  small  parameter  (electron/ion  mass  ratio)  in  the  plasma 
was  outlined.  We  prepared  detailed  work  plan  for  application  of  these  methods  during 
Phase  II  research. 

•  We  have  tested  different  Eulerian  and  semi-Lagrangian  kinetic  solvers  for  collisionless 
plasmas.  An  efficient  high-order  Vlasov-Poisson  solver  has  been  demonstrated  for  self- 
consistent  simulation  of  a  collisionless  boundary  layer.  Better  understanding  of  the 
collisionless  particle  kinetics  in  a  DC  sheath  was  achieved  as  reflected  in  the  book 
chapter  submitted  for  publication. 

•  We  have  developed  a  local  Boltzmann  solver  with  adaptive  mesh  in  velocity  space  and 
tested  this  solver  for  simple  collision  integrals  describing  electron  and  ion  collisions  in 
weakly  ionized  plasmas. 

•  We  have  developed  hydrodynamic  models  of  plasmas  with  octree  based  dynamically 
adaptive  Cartesian  mesh  and  demonstrated  the  model  capabilities  for  streamer 
simulations. 

•  The  paper  "Streamer  Simulations  with  Dynamically  Adaptive  Cartesian  Mesh"  was 
prepared  and  accepted  for  publication  in  a  special  issue  of  IEEE  Transactions  on 
Plasma  Science  devoted  to  "Images  in  Plasma  Physics". 
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1. 


APPLICATION  OF  INVARIANT  MANIFOLD  AND  RENORMALIZATION 
GROUP  METHODS  TO  PLASMA  MODELING 


The  renormalization  group  (RG)  method  has  been  recently  demonstrated  as  a  powerful  tool  for 
reduction  of  evolution  equations  in  terms  of  invariant  manifolds.  It  was  recognized  [1]  that  the 
reduction  of  evolution  equations  is  a  natural  extension  of  the  well-known  asymptotic  method 
by  Krylov,  Bogoliubov  and  Mitropolski  for  nonlinear  oscillator  -  the  later  is  in  fact  an  RG 
theory  although  the  term  RG  is  not  used.  Using  this  methodology,  the  Boltzmann  equation  was 
derived  as  an  RG  equation  [2],  the  Navier-Stokes  equations  were  derived  from  the  Boltzmann 
equation  [3],  and  stable  post-Navier-Stokes  equations  have  been  obtained  [4]  instead  of  the 
unstable  Burnett  equations.  For  a  classical  viscous  flow  problem,  the  RG  calculations  have 
lead  to  a  new  prediction  for  the  drag  coefficient,  which  can  reproduce  and  surpass  the  results  of 
matched  asymptotics  [5], 

1.1  Slow  Invariant  Manifolds 


The  problem  of  model  reduction  in  physical  kinetics  was  recognized  as  a  problem  of  time 
separation  and  construction  of  slow  invariant  manifolds.  To  illustrate  these  notations,  consider 
an  //-dimensional  dynamical  system  governed  by  the  evolution  equation 

<iW 

—  =  F(W,0  (1) 

dt 

When  the  dynamics  is  reduced  to  an  m-dimensional  system  ( m<n ),  the  vector  W(t)  approaches 
a  well-defined  m-dimensional  manifold  M  embedded  in  the  ^-dimensional  phase  space  (see 
Figure  1).  The  geometrical  object  Mis  called  an  attractive  manifold.  If  after  some  time  W(t)  is 
confined  in  the  manifold  M,  then  M  is  called  an  invariant  manifold.  Furthermore,  if  the 
dynamics  on  M  is  slow,  M  is  called  a  slow  manifold.  There  is  no  explicit  definition  of  slow 
invariant  manifolds  without  explicit  small  parameter  in  the  system.  It  was  demonstrated  how 
RG  method  helps  construct  invariant  manifolds  and  give  reduced  dynamics  on  them  [3].  The 
concept  of  invariant  manifolds  for  physical  and  chemical  kinetics  was  explicitly  developed 
very  recently  [6,7]. 


/  *  -  to  ru 

V 


M  =  m  <  n 
” - - 

Figure  1.  Illustration  of  the  reduced  dynamics.  The  dynamic  variable  W(t)  approaches  to 
and  after  some  time  is  eventually  confined  in  the  manifold  M  of  lower  dimensionality. 

1.2  Collisionless  Plasma 

During  Phase  I,  we  have  studied  how  the  RG  and  invariant  manifold  concepts  can  be  applied  to 
modeling  weakly  ionized  plasmas.  We  first  tried  to  apply  strict  results  of  quasi-neutral  plasma 
theory  to  interpret  dynamics  of  the  plasma-wall  transition  layer  (presheath  in  CCP,  and  skin 
layer  in  ICP).  The  analytical  expression  for  the  spatio-temporal  evolution  of  the  distribution 
function  of  electron  and  ions  having  initial  Gaussian  distribution  in  space  and  Maxwellian 
distribution  in  velocities  has  been  obtained  by  RG  methods  in  the  form 
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(2a) 


P  = 


PPn\i. 


f  = 


-exp 


-exp 


1  +  r2 

\v-uf 

(l  +  m/ M)  eol 

VTe 

2 

V 

l-Vxi /<  m  J 

l  +  T 


(v-m)2  (l  +  //2)  eo 
2  1  —  VTe^VTi  171 


(2b) 


where  the  local  mean  velocity  of  the  plasma  u  and  the  electric  potential  ®  are  defined  as 


l  +  T2 


®: 


2e(l  +  r2)2  (1  +  p2) 


Here  r  =  Clt,  g  =  x/L,  p  =  -> JmJM  ,  cs  =  ^j(Ti  +T,)/(M  +  m)  is  the  sound  velocity  based  on  the 
initial  temperatures  of  electron  and  ions,  Q  =  cs  /  L,  and  Lisa  characteristic  scale  of  plasma 
non-uniformity  (not  necessarily  equal  to  the  Debye  radius  rDe).  This  solution  is  valid  for 
arbitrary  ratio  of  the  electron  and  ion  mass,  ml  M  ,  and  arbitrary  values  of  the  temperatures,  Te 
and  T. . 


The  analytical  solution  (2)  describes  many  features  observed  in  our  simulations  of  the 
collisionless  plasma-wall  transition  layer  (see  below)  and  results  of  other  authors.  In  particular, 
during  the  time  evolution  of  the  initial  perturbation,  a  shift  of  the  mean  velocity  of  electron  and 
ions  takes  place,  and  the  temperature  of  ions  decreases.  Due  to  the  difference  of  electron  and 
ion  mass,  M  >m  ,  ion  acceleration  is  more  pronounced  since  for  the  thermal  velocities  vTi<vTe, 
and  the  electron  acceleration  is  negligible. 

During  Phase  I,  we  obtained  approximate  solutions  of  unsteady  Vlasov-Maxwell  system  using 
RG  algorithms  and  slow  invariant  manifolds.  We  considered  one-dimensional  system 
described  by  the  system  of  kinetic  equations 


d+v 

dt  “  dx  ma  dva 


=  0 


a  =  e,i 


(3) 


with  additional  equations  for  the  electric  field 


dE  .  .  A  dE  . 

—  +  4xj=0,  - — 4np  =  0, 

at  ox 


(4) 


and  nonlocal  material  relations 

5?  +  £  =  0’  P  =  Hea\fad^a’  j  =  Y,ea\fPad^a  (5) 

We  obtained  solutions  of  this  system  for  several  cases  described  below. 

Example  1.  This  case  corresponds  to  generalization  of  the  quasineutral  solution  (2)  for  non¬ 
zero  initial  velocity  of  the  plasma  constituents  having  Maxwellian  velocity  distribution 
functions.  Such  a  solution  can  be  obtained  by  applying  Galilean  transformation  to  the  solution 
obtained  previously.  From  the  physical  stand  point  this  operation  is  equivalent  to 
transformation  to  a  moving  reference  frame,  obtaining  solution  in  this  frame,  and  then 
returning  to  the  laboratory  reference  frame.  The  velocity  distribution  functions  and  electric 
field  in  the  plasma  are  given  by 
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exp 


(6a) 


r 


\f2aVj 


1  +  TZ 


2VT  a  V 


cX  +  U  Y  -UT  /  Cs)2 


1  +  r 


2(l  +  r2) 


fc-«r/c,)(v*e-v*e),  (6b) 

(l  +  r2)  (l  +  w^/m,.) 

It  follows  from  this  solution  that  the  mean  particle  velocity  at  x=0,  which  is  equal  to  u  at  t=0, 
decrease  monotonically  with  increasing  time  t.  The  particle  density  reaches  maximum  value  at 
the  pointy  =  ut/cs,  which  is  moving  with  time.  It  is  exactly  at  this  point  that  mean  velocity 

keeps  its  initial  value,  u.  In  fact,  the  solution  (6)  describes  a  localized  perturbation  moving 
across  the  plasma  with  velocity  u. 

The  solution  (6)  can  be  obtained  by  other  means.  An  alternative  approach  consists  in  seeking  a 
solution  to  kinetic  equations  that  is  invariant  with  respect  to  a  group  of  local  transformations 
with  a  generator,  which  appears  as  a  linear  combination  of  generators  of  space  and  time 
translations  and  the  projective  group  generator. 

The  solution  (6)  is  a  particular  case  of  a  more  general  solution,  which  does  not  assume 
Maxwellian  VDF.  In  the  general  case,  the  dependence  of  the  electric  field  on  spatial  coordinate 
may  deviate  from  the  linear  dependence  (6b).  This  general  solution  has  the  form 


fa 


r 


cX+«t  ,  {c£~UTf  , 


1  +  r 


+ 


2(l +  r2)2 


+  -^0(7) 


(7a) 


E 


1  d<S> 

(1  +  r2)3/2  U’ 


{c£-ur) 

VT77  ’ 


(7b) 


The  dependence  of  O(J)  is  determined  by  the  initial  VDFs  of  the  particles  and  by  the 

condition  of  plasma  quasineutrality.  The  quantities  Ia  and  J  in  (7)  are  invariants  of  the  group 
generator  constructed  from  the  generators  of  time  and  space  translations  and  the  projective 
group  generator.  The  VDFs  do  not  change  under  group  transformation,  i.e.  they  remain 
invariant,  and  their  relation  to  the  initial  VDFs  expressed  through  Ia  and  J  defines  an 
invariant  manifold  corresponding  to  the  solution  of  the  problem. 

Example  2.  As  shown  in  the  previous  example,  the  group  of  Galilean  transformations 
transforms  any  solution  of  the  kinetic  equations  with  zero  mean  velocity  to  a  solution  with 
non-zero  mean  velocity.  As  stated  above,  one  can  obtain  the  same  solution  by  an  alternative 
approach  consists  in  seeking  a  solution  to  kinetic  equations  invariant  under  the  local 
transformation  group  with  a  generator,  which  is  a  linear  combination  of  generators  of  time  and 
space  translations;  in  the  general  case  one  can  also  add  the  generator  of  the  Galilean  group. 
Since  all  these  group  generators  are  admitted  by  the  original  Vlasov-Maxwell  system  without 
quasineutrality  assumption  used  for  obtaining  (6),  the  solution  (8)  below  is  valid  for  arbitrary 
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relation  between  the  Debye  length  and  the  characteristic  spatial  scale  of  the  plasma 


fa  =Fa{Ca)  Ca  =-(v0 -u-at)2  +aA+  -^O(A), 
2 v  ’ 


(8a) 


E  =  — 


JO 

dA’ 


A  =  x  —  ut  — 


at 


(8b) 


The  dependence  of  O(A)  is  found  from  the  initial  VDFs  with  the  help  of  Poisson  equation. 
The  solution  (8)  is  a  stationary  solution  of  the  Vlasov  kinetic  equations  obtained  in  a  reference 
frame,  which  moves  with  initial  velocity  u  and  acceleration  a.  The  solution  (8)  has  similar 
group  interpretation  as  the  solution  (7),  namely  C“  and  A  are  invariants  of  the  corresponding 
group  generator,  which  describes  the  solution  to  the  problem. 

The  obtained  solutions  to  the  one-dimensional  Vlasov-Maxwell  system  are  quite  obvious  from 
the  physical  standpoint  and  thus  could  be  obtained  earlier.  However,  specific  solutions  to 
particular  problems  could  be  very  useful  from  different  points  of  view.  In  particular,  the 
developed  procedure  can  be  applied  to  two-component  plasmas  with  prescribed  ion  motion. 
The  ions  can  be  assumed  immobile  with  given  distribution  of  density,  or  moving  in  space  with 
a  given  velocity.  The  corresponding  expressions  can  be  obtained  from  the  general  formulas. 

It  appears  to  be  difficult  to  obtain  general  solutions  for  arbitrary  initial/boundary  conditions. 
However,  using  the  concept  of  slow  invariant  manifolds,  one  can  expect  to  construct  an 
approximate  symmetry  (RG  type)  for  the  Vlasov-Maxwell  system,  which  would  allow  one  to 
construct  solutions  for  quite  arbitrary  initial  conditions.  However,  such  a  procedure  would 
require  serious  research  work. 

1.3  Collisional  Plasma 

During  Phase  I,  we  have  also  studied  how  invariant  manifold  methods  could  be  used  for 
reduced  description  of  weakly  ionized  plasmas  with  collisions  among  particles.  The  difference 
between  scattering  rate  and  the  energy  loss  rate  for  electrons  moving  in  a  neutral  gas  under 
effect  of  external  electric  field  could  be  utilized  for  the  reduced  description  of  electron  kinetics. 
For  slow  electrons  with  energies  less  than  the  excitation  potential  of  atoms  (~10  eV),  scattering 
dominates  over  energy  loss.  For  these  electrons,  the  reduced  description  is  possible  along  the 
following  lines:  Boltzmann  kinetic  equation  — >  two-term  Spherical  Harmonics  Expansion 
(Fokker-Planck  Equation),  nonlocal  approach  — »  hydrodynamic  model.  For  fast  electrons,  the 
relative  importance  of  scattering  and  energy  loss  can  be  opposite  compared  to  slow  electrons. 
For  fast  electrons,  small  angle  scattering  dominates  resulting  in  electron  beams  and  runaway 
electrons  in  high  electric  field.  We  anticipate  that  the  Invariant  Manifolds  and  RG  methods  can 
be  applied  for  the  renormalization  of  the  Boltzmann  collision  integral  and  reduced  description 
of  collisions  in  terms  of  nonlocal  friction  force  that  depends  on  the  distribution  function. 

We  have  identified  several  small  parameters  that  can  be  used  for  reduced  description  of  the 
system.  These  are 

•  the  ratio  of  electron  to  ion(atom)  mass  me  /  ma 

•  the  ratio  of  plasma  to  gas  density,  i]  =  ne  /  N  (ionization  degree) 

•  the  ratio  of  inelastic  to  elastic  collision  frequencies  for  electrons,  v  !v . 
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The  first  one  is  probably  the  most  important  and  defines  key  features  of  gas  discharge  physics. 
As  a  result  of  the  large  difference  between  the  electron  and  atom  mass,  electrons  respond 
nearly  adiabatically  to  changes  of  the  external  electromagnetic  fields,  and  plasma  shields  the 
fields  with  the  (fast)  electron  time  scale.  The  ion  transport  and  ionization  processes  occur  with 
(slow)  ion  time  scale.  Collisionless  invariants  (such  as  total  energy  (kinetic  plus  potential))  are 
important  characteristics  of  the  Vlasov-Maxwell  system.  Small  deviations  from  the  adiabatic 
electron  motion  can  be  responsible  for  the  collisionless  (stochastic)  electron  heating  -  the 
concept  proposed  for  detailed  Phase  II  studies. 

Small  ratio  of  inelastic  to  elastic  collision  frequencies  for  electrons  at  energies  of  the  order  of 
the  excitation  potential  of  atoms  (~  1 0  eV)  is  another  important  parameter  enabling  reduced 
description  of  electron  kinetics  in  gas  discharges.  As  a  result  of  me/  maU  1  and  v*  /  v  □  1 ,  the 

reduced  description  of  electrons  can  proceed  along  the  following  lines:  Boltzmann  Transport 
Equation  ->  two-term  Spherical  Harmonics  Expansion  (Fokker-Planck  Equation),  nonlocal 
approach  ->  hydrodynamic  model  (see  Figure  2).  The  Fokker-Planck  approach  to  simulation  of 
electron  kinetics  in  collisional  gas  discharge  plasmas  was  implemented  in  the  commercial 
software  for  plasma  simulations  developed  at  CFDRC  [8].  Many  aspects  of  the  kinetic  theory 
of  gas  discharges  [9]  can  be  interpreted  in  terms  of  invariant  manifolds.  In  particular,  the 
procedure  of  spatial  averaging  of  Bemstein-Holstein-Tsendin  for  calculation  of  the  electron 
distribution  function  in  DC  and  RF  discharges  (nonlocal  approach)  can  be  recognized  as  a 
constructive  method  for  creating  slow  invariant  manifold  for  electron  kinetics. 


Figure  2.  A  hierarchy  of  transport  models  for  electrons  (from  [8]) 


For  fast  electrons,  the  relative  importance  of  scattering  and  energy  loss  can  be  quite  opposite 
compared  to  slow  electrons.  Possible  scenarios  of  EDF  formation  in  a  steady  electric  field  E 
for  arbitrary  rates  of  scattering  and  deceleration  depicted  qualitatively  in  Figure  3  taken  from 
Ref.  10. 
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In  the  domain  a  both  the  EDF  body  and  its  tail  are 
almost  isotropic.  The  random  walk  in  energy  occurs 
with  a  step  eEA ,  which  is  small  with  respect  to  the  EDF 
characteristic  scale  e1,  so  one  can  use  the  concept  of 
diffusion  in  energy  with  a  diffusion  coefficient 
D..(v)  =  (eEX)2  v/3 .  In  the  domain  b  the  energy  loss 

can  be  treated  as  quasi-elastic  and  the  distinction 
between  the  EDF  body  and  tail  disappears.  In  the 
domain  c  the  random  walk  step  eEA  is  small,  therefore, 
the  EDF  body  can  be  approximated  by  an  isotropic 
"pipe-line"  EDF.  In  the  tail,  however,  inelastic  collisions 
are  very  frequent,  resulting  in  a  large  anisotropy  with  no 
electrons  moving  against  the  field.  In  the  domain  d,  both 
the  EDF  body  and  tail  are  strongly  anisotropic,  and  a 
needle-like  EDF  along  the  electric  field  is  formed  at  energies  w  >  sx .  A  small  isotropic  halo 

from  elastically  scattered  electrons  contains  a  small  fraction  of  the  total  number  of  electrons. 
Fast  electrons  experience  small  angle  scattering  and  can  penetrate  deeply  into  a  dense  gas 
producing  non-local  ionization  in  the  areas  with  no  electric  field.  The  problem  of  fast  electron 
kinetics  in  atmospheric  pressure  plasmas  has  recently  attracted  increased  attention  due  to 
numerous  technological  applications.  We  have  previously  developed  simplified  models  for  the 
fast  electron  transport  and  have  identified  the  problem  of  runaway  electrons  as  a  target  area  for 
Phase  II  research.  It  is  anticipated  that  the  Invariant  Manifolds  and  RG  methods  can  also  be 
applied  for  the  renormalization  of  the  Boltzmann  collision  integral  and  reduced  description  of 
collisions  in  terms  of  nonlocal  friction  force  that  depends  on  the  distribution  function.  An 
example  of  such  a  description  can  be  found  in  [1 1], 


Figure  3  Electron  kinetics  in 
external  electric  field  for  arbitrary 
ratios  of  scattering  and 
deceleration 
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2. 


KINETIC  SOLVERS  FOR  COLLISIONLESS  PLASMAS 


We  have  tested  accuracy  and  performance  of  several  Eulerian  and  semi-Lagrangian  kinetic 
solvers.  We  have  applied  a  high-order  Vlasov-Poisson  solver  for  self-consistent  simulation  of  a 
boundary  layer  in  collisionless  plasma. 

2.1  Eulerian  Kinetic  Solvers 

Three  Vlasov  solvers  have  been  tested  for  a  simple  case  related  to  DC  sheath  in  gas  discharge 
plasmas.  The  first  solver  is  an  extension  of  the  UFS  kinetic  solver  [12]  with  added  electric 
force.  The  second  solver  was  generously  provided  by  Dr.  Shoucri  [13].  The  third  solver  is  the 
Convective  Scheme  solver  described  in  [14].  We  considered  electric  field  in  the  form 


E(x)  = 


iA(  sm~x)/sm, 

I  o 


0  <  x  <  sm 


X  >  s 


m 


(9) 


which  corresponds  to  the  potential  cp{x)  =  -Asm(\-xl sm)\\-(\  +  xl sm)!2\ .  The  particles  are 

injected  at  x  =  sm  with  a  velocity  distribution  function 

f(vx ,  vv )  =  exp  |- ( vx  -  E(l  )2  /  2 j  exp  (-vv2 12)  and  either  accelerated  or  decelerated  by  the 

electric  field  in  the  sheath  (depending  on  sign  of  the  electric  charge).  The  analytical  solution 
for  the  VDF  has  been  found  in  the  form 


f(vx,vv)  =  exp 


{I 


'  2(p{x)  -  V0 


2  ^ 

12 

J 


exp(-v//2) 


(10) 


According  to  this  solution,  the  ion  VDF  in  the  plane  (vx,vv)  is  deformed  (become  narrower 
along  vx)  during  ion  acceleration  by  the  DC  field.  This  means  that  the  ion  temperature  T 
decreases  in  the  sheath. 


Figure  4.  UFS  ion  density  and  temperature  from  (on  the  left)  for  different  size  of  velocity 

and  space  mesh.  Ion  temperature  from  ShoucrVs  code  and  UFS  (on  the  right) 
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Figure  4  (left  part)  shows  the  UFS  numerical  solution  for  =  90 ,  F0  =  10 ,  A  =  8 .  This  force 

corresponds  to  the  potential  (p(Q)  =  360  and  the  mean  velocity  is  v(0)  =  V820  =  28.6 .  It  is  seen 
that  even  with  the  spatial  mesh  of  200  points  and  velocity  mesh  of  250  points  we  can  not 
reproduce  the  analytical  profile  of  the  ion  temperature  with  the  UFS  code.  The  right  part  of 
Figure  4  shows  the  temperature  profiles  obtained  with  the  UFS  and  S-codes.  The  UFS  code 
uses  second  order  scheme  with  van  Leer  limiter,  the  S-code  uses  cubic  interpolation.  The  S- 
code  reproduces  analytical  result  even  for  (100,120)  nodes.  The  UFS  fails  even  for  (200,250) 
nodes. 

Figure  5  shows  velocity  distributions  from  S-code  and  UFS  for  potential  drop  400.  The  UFS 
results  are  for  200  nodes  in  physical  space  and  250  nodes  in  velocity  space.  The  S-results  are 
for  (200,  240),  (100,120)  and  (50,60)  nodes. 


Figure  5.  Ion  distribution  functions  from  Shoucri's  code  (left)  and  UFS  (right). 

For  V0  =  0 ,  the  velocity  distribution  remains  isotropic,  and  the  temperature  should  be  constant 

in  the  sheath.  Figure  6  shows  the  electron  distribution  function  at  the  electrode  and  temperature 
distributions  from  the  Shoucri's  and  UFS  codes.  The  potential  drop  is  -8.  The  UFS  results  are 
for  200  nodes  in  physical  space  and  100  nodes  in  velocity  space.  The  S-results  are  for  50  nodes 
in  space  and  140  nodes  in  velocity  space.  It  is  seen  that  the  S-code  reproduces  the  expected 
temperature  profile  with  better  accuracy.  The  temperature  drop  near  the  boundary  is  due  to 
absorption  of  the  fastest  electrons  at  the  boundary. 
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Figure  6.  Electron  distribution  functions  from  Shoucri's  code  at  electrode  (left)  and 
spatial  distribution  of  electron  temperature  (right). 


2.2  DC  boundary  layer  in  collisionless  plasma 

We  have  derived  analytical  expressions  for  the  electron  and  ion  velocity  distribution  functions 
in  the  collisionless  DC  sheath  and  verified  results  of  the  numerical  solution  of  the  Vlasov- 
Poisson  system  [15]  versus  the  analytical  solution. 

The  problem  statement  is  described  in  details  in  [15].  Briefly,  electrons  and  ions  are  injected 
from  the  right  boundary  (x=L)  with  given  velocity  distributions.  The  boundary  at  x=0 
corresponds  to  floating  potential  collecting  charge.  We  use  dimensionless  variables  introduced 
in  [15],  where  space  is  normalized  to  the  electron  Debye  length,  time  is  in  units  of  inverse  ion 

plasma  frequency,  and  velocity  is  in  the  units  of  ion  acoustic  speed,  cs  =  ^Te  /  M  .  A  surface 

charge  is  built  up  at  x=0  to  equalize  electron  and  ion  currents  to  the  wall  in  a  steady  state.  The 
problem  is  characterized  by  two  parameters:  the  ratio  of  electron  to  ion  mass,  xmei=m/M,  and 
the  ratio  of  electron  to  ion  temperatures,  tei  =TeITr 

In  our  present  simulations,  L- 20,  ion  time  step  At  =  0.01.  We  solve  electron  Vlasov  equation 
with  Poisson  equation  in  a  sub-cycle  with  nturn  steps  and  the  electron  time  step  At  /  ntiirn 
(with  nturn= 60).  The  velocity  space  for  electrons  -4 /  yjxm~i  <  d,e  <  4 /  A Jxm ~ei  ,  for  ions 

-17/ yjt~  <  f  <  0 / yft~  at  tei  =  30.  The  problem  reaches  a  steady  state  at  ?~10  in  less  than  an 
hour  of  computing  time.  Results  shown  below  are  for  t= 20. 

Figure  7  shows  spatial  distributions  of  particle  densities  and  temperatures  for  tei= 30.  The 
electron  density  drops  in  the  sheath  because  the  electrostatic  potential  repels  the  electrons  back 
to  plasma.  The  analytical  solution  of  the  Vlasov  equation  for  the  electron  velocity  distribution 
function  has  the  form 


/«(£*) 


exp(-£')  g  <  ^2A(p(x)  /  xmei 

0  £>72A^(x)/  Xmei 


(11) 
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where  s  - 1;1  /  2  -  cp(x)  is  the  total  energy,  cp(x)  is  the  electric  potential  measured  in  units  of 
Te  ,  and  Acp(x)  =  (p{G)-(p{x)  .  It  is  assumed  that  (p(L)  =  0.  The  electron  temperature  is 

dropped  in  a  vicinity  of  the  wall  due  to  the  depletion  of  electron  velocity  distribution  by  fast 
electrons,  which  overcome  the  electrostatic  potential  barrier  and  escape  to  the  wall.  No 
electrons  reaching  the  wall  return  back  to  plasma,  and  the  electron  velocity  distribution 

function  is  zero  at  ^  >  ^2A cp(x)l  xmei  (see  Figure  8). 


The  ion  density  drops  because  the  electrostatic  field  accelerates  the  ions  and  their  mean 
velocity  increases  while  the  ion  current  remains  constant.  The  ion  velocity  distribution  can  be 
described  by  the  analytical  expression 

f  ,  , -  \2  A 


exp 


-t 


■i(\[¥+  2p(*)-v„)  / 2 


£<0 

£>0 


(12) 


It  follows  from  this  expression  that  the  amplitude  of  the  ion  distribution  remains  constant,  and 
the  width  of  with  respect  to  velocity  decreases  in  the  sheath.  The  decrease  of  the 


IVDF  width  corresponds  to  decrease  of  the  ion  density  and  temperature  Tix(x)  in  the  sheath 
(see  Figure  7). 


Figure  7.  Particle  densities  (left)  and  temperatures  (right)  for  tei=30. 
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Figure  8.  Electron  velocity  distributions  (left)  and  ion  velocity  distributions  (right)  at 
different  spatial  positions  for  tei=30.  Numbers  denote  spatial  mesh  points  for  the  total  of  400 

mesh  points. 


Figure  9  compares  the  distribution  of  density  profiles  and  temperature  profiles  for  tei=\0  and 
tei= 30.  It  is  seen  that  the  sheath  width  is  slightly  larger  at  tei=  10  compared  to  tei= 30. 


Figure  9.  Comparison  of  density  profiles  (left)  and  temperature  profiles  (right)  for  tei=10  and 

tei=30. 


Figure  10  shows  the  distribution  functions  for  Ar,  which  can  be  compared  with  Figure  8  to 
understand  changing  of  the  velocity  range  with  changing  xmei  parameter. 
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Figure  10.  Comparison  of  density  profiles  (left)  and  temperature  profiles  (right)  for  tei=10  and 

tei=30. 


Our  research  described  in  this  section  allowed  better  understanding  of  particle  kinetics  in  the 
collisionless  DC  sheath.  The  results  were  reflected  in  a  book  chapter  [16]  submitted  for 
publication. 

2.3  Convected  Scheme 

During  Phase  I,  the  academic  partner  at  the  University  of  Wisconsin  has  tested  the  Convected 
Scheme  (CS)  method  for  collisionless  sheath  problems.  The  CS  is  a  semi-Lagrangian  solution 
to  the  Boltzmann  equation,  which  offers  advantages  over  standard  methods  of  characteristics  in 
that  it  naturally  conserves  density,  energy  and  other  variables,  locally  in  phase  space.  It  is 
believed  to  be  more  accurate  than  fully  mesh-based  (Eulerian)  methods.  The  CS  offers 
advantages  over  other  methods  in  the  "intermediate”  regime  where  the  mean  free  path  is 
comparable  to  the  spatial  scales,  due  to  its  ability  to,  on  one  hand,  take  steps  up  to  a  fraction  of 
the  mean  free  path,  and  on  the  other,  provide  conservation  at  each  step  it  takes. 

As  a  first  step  towards  hybrid  models,  we  have  set  up  a  self-consistent  and  time-dependent 
model  for  the  electron  and  ion  distribution  functions  in  the  RF  sheath.  In  Figure  11  we 
compare  CS  results  for  the  electron  distribution  to  an  analytical  solution  for  a  given 
distribution  of  the  electrostatic  potential;  very  close  agreement  is  obtained.  The  results  are 
consistent  with  Figure  8  obtained  with  the  Eulerian  Vlasov  solver,  which  is  especially 
revealing  when  the  velocity  distributions  are  plotted  versus  electron  velocity  rather  than  kinetic 
energy  (see  Figure  11). 
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Figure  11  Electron  distribution  functions  plotted  versus  kinetic  energy  (left)  and  velocity 

(right)  for  different  points  inside  the  sheath. 


The  time  evolution  of  the  ion  and  electron  densities  in  the  collisionless  sheath  at  RF  frequency 
200  MHz  is  shown  in  Figure  12.  In  these  calculations,  the  electron  density  follows  the 
Boltzmann  law,  and  the  Vlasov  equation  for  ions  is  solved  together  with  Poisson  equation  with 
boundary  conditions  for  the  potential  <p{Qi)  =  %  (1  +  cos(&tf ))  *=0  and  cp(L)  =  0  .  cpQ  =  200  V. 


Number  of  ion  (solid)  and  electrons  (dashed)  in  the  sheath 
Sheath  width  -  0.375  turn 


Figure  12  Time  evolution  of  the  ion  and  electron  densities  in  the  sheath 


In  the  future,  we  intend  to  run  kinetic  calculations  for  both  the  ion  and  electron  distributions  in 
the  RF  sheath,  and  develop  methods  to  accelerate  simulations  taking  advantage  of  the  disparity 
of  electron  and  ion  mass.  We  will  then  perform  analysis  of  the  two  sets  of  results  obtained,  to 
ensure  that  the  CS  solution  for  the  electron  distribution  matches  that  obtained  using  the 
simplified  relation. 
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3. 


BOLTZMANN  SOLVER  WITH  ADAPTIVE  MESH  IN  VELOCITY  SPACE 


It  is  known  that  Lagrangian  methods  are  more  efficient  than  Eulerian  methods  in  terms  of 
computer  memory  because  they  use  only  those  parts  of  phase  space  where  the  particles  are 
present.  The  efficiency  of  Eulerian  kinetic  solvers  can  be  increased  by  using  phase  functions  or 
adaptive  mesh  in  velocity  space.  In  fact,  some  people  believe  that  due  to  high  dimension  of 
phase  space,  adaptive  mesh  is  a  must  for  multi-dimensional  kinetic  simulations  using  Euler 
methods.  Octree  Cartesian  mesh  is  particularly  attractive  for  this  purpose  because  there  are  no 
complex  boundaries  in  velocity  space. 


We  have  implemented  local  Boltzmann  solver  with  adaptive  mesh  in  velocity  space  for  simple 
types  collision  terms.  The  local  kinetic  equation  contains  two  parts:  the  transport  in  velocity 
space  under  action  of  an  external  force  and  the  collisional  part: 


cf  [  eE(t)  df 
dt 


m 


34 


=5 


(13) 


The  transport  part  was  implemented  (in  2D  and  3D  velocity  space)  using  the  finite  volume 
formulation  with  second-order  accuracy  in  time  and  space.  The  octree/quadtree  data  structure 
allowed  dynamic  grid  adaptation.  For  mapping  from  one  velocity  grid  to  another,  second-order 
accuracy  scheme  was  used.  The  collisional  part  depends  on  the  type  of  the  collisional  integral. 
We  have  implemented  so  far  three  types  of  collision  integrals 


•  the  BGK  collisional  integral 

•  the  charge-exchange  collisional  integral  for  ions  moving  in  a  parent  gas 

•  the  elastic  (isotropization)  part  of  the  collision  integral  for  electrons  collisions  with 
atoms  (in  the  limit  mIM  — >  0). 

For  the  BGK  collisional  integral,  the  use  of  adaptive  meshes  allowed  us  to  avoid  special 
corrections,  which  required  before  for  coarse,  non-adaptive  grids.  The  isotropization  integral 
for  electrons  was  implemented  using  a  special  technique  developed  in  our  previous  work  to 
correctly  treat  leaps  of  particles  from  one  phase  cell  to  another  on  non-uniform  and 
unstructured  grids. 

All  calculations  presented  below  were  very  fast,  taking  from  30  sec  to  10-20  min  on  a  2.2  GHz 
PC.  This  is  due  to  the  adaptive  grid  capabilities  allowing  us  to  minimize  the  number  of  velocity 
grid  points  to  adequately  describe  the  studied  distribution  functions. 

3.1  BGK  Collision  Integral 


We  first  tested  a  collisionless  case  for  a  time  varying  electric  field,  E  =  Eocos(cot).  We  have 
confirmed  that  the  density  and  temperature  are  conserved  with  high  precision  during  the 
simulations,  while  the  mean  velocity  oscillated  with  the  field  frequency  to.  Then,  the  effect  of 
collisions  was  investigated.  The  BGK  collision  term  was  used  in  the  form  v(fm  -  f)  where  fm 

is  a  Maxwellian  distribution  with  prescribed  density,  temperature  and  mean  velocity,  and  v  is 
the  collision  frequency.  Figure  13  shows  the  computational  mesh  and  2D  distribution  function 
in  velocity  space  (4x,4y) .  Once  can  see  that  the  velocity  grid  adapts  on  the  features  of  the 

distribution  function  based  on  gradient  of  the  velocity  distribution  function.  Figure  14  shows 
the  velocity  distributions  along  2,  for  2,  — 0  for  two  different  times.  The  final  distribution 

function  (after  about  500  collisional  times)  has  two  peaks.  These  peaks  correspond  to  positions 
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in  velocity  space  where  the  particles  spend  most  of  the  time  (where  dE/dt  is  small  and 
changing  sign). 


Figure  14.  The  velocity  distribution  function  with  respect  of  for  f  =  0 . 

Using  adaptive  mesh  in  velocity  space  allows  accurate  description  of  the  velocity  distribution 
function  and  its  moments.  An  analytic  expression  for  the  mean  velocity  was  obtained  in  a 
simple  form: 

u  =  u0  exp(-U)  +  um  [1  -  exp(-vt)]  +  E0  ![v2  +  of  ]  [v  cos  (cot)  +  co  sin(<yt)  -  v  exp(- vt)] 


where  u0  is  the  initial  velocity.  The  equation  for  temperature  was  solved  numerically. 

The  comparison  of  moments  calculated  from  VDF  with  analytical  results  demonstrated  very 
good  agreement,  thus  proving  that  the  numerical  scheme  is  capable  of  maintaining  high 
accuracy  during  simulations  of  long  time  evolution  (here,  about  500  collisional  times).  In 
addition,  we  have  checked  that  the  implemented  BGK  collision  integral  maintains  the  shape  of 
the  distribution  function  in  the  c,  -direction.  The  test  showed  that  the  Tv  temperature  remains 

constant  over  the  simulated  time  interval  with  a  2%  precision. 

3.2  Ions  in  a  DC  electric  field  with  charge-exchange  collisions 

For  charge-exchange  collisions  of  ions  with  a  parent  gas  the  collision  integral  can  be  written 
in  the  form  [17] 

5,.  =  n  J  1 5-5, 1  ^1  (14) 

R 3 

where  <p(  f)  is  a  Maxwellian  velocity  distribution  of  the  background  gas  atoms  with  density  N 
and  temperature  T,  crex  is  a  charge  exchange  collision  cross-section. 
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Consider  ions  with  a  mass  M  injected  in  DC  electric  field  E  and  experiencing  charge-exchange 
collisions  with  a  steady  background  gas.  The  charge  exchange  collisions  cannot  change  the 
velocity  distribution  orthogonal  to  the  electric  field  direction.  So,  the  velocity  distribution  is  of 
the  form  /(^)  =  <f>(4r1)F(4rx) .  The  distribution  function  F(£x)  evolves  in  time  until  a  steady 


state  is  reached.  The  analytical  expression  for  the  steady  state  distribution  F(£x)  can  be 
obtained  in  the  form 


F(£)  =  Cex p 


NM 

~eE 


(15) 


Figure  15  shows  results  of  the  numerical  solution  of  the  2D  kinetic  equation  with  adaptive 
mesh  in  velocity  space  for  the  ion  distribution  function  under  an  applied  force  of  magnitude  4 
and  charge-exchange  collision- frequency  of  1.  Figure  15  shows  the  calculated  ion  distribution 
function  at  three  moments  in  time:  initial,  intermediate  and  final.  One  can  see  the  grid 
adaptation  on  the  features  of  the  distribution  function  as  it  evolves  in  phase  space.  Figure  16 
shows  the  distribution  functions  along  the  =  0  line  for  different  time  moments  (left)  and 

distribution  functions  along  =const  lines.  One  can  see  that  the  calculated  steady-state 


distribution  function  compares  well  with  the  analytical  solution.  In  the  £  -direction,  the 

distribution  function  remains  a  Maxwellian  with  the  temperature  of  the  source  function  (which 
represents  cold,  motionless  atoms  on  which  charge-exchange  occurs). 


final. 


Figure  16.  Ions  distribution  functions  along  d  =  0  line  for  different  time  moments  (left)  and 
distribution  functions  vs  d  lines  for  different  %x  positions  (right). 


Figure  17  shows  evolution  of  macroparameters  of  the  simulated  ion  distribution  function  in  an 
external  field  under  charge-exchange  collisions.  One  can  see  that  the  density  is  well  conserved 
(variation  is  less  than  0.05%),  the  velocity  and  temperatures  reach  steady  state  in  around  5-6 
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collision  times.  One  can  also  notice  that  the  temperature  (or  average  energy)  in  the  field 
direction  (here,  Tx)  becomes  about  300  times  larger  than  the  transversal  one  (here,  Ty ),  which 
itself  remains  constant  (recall  that  the  charge-exchange  collision  integral  does  not  change  the 
distribution  function  in  the  transversal  direction).  These  facts  provide  evidence  of  that  the 
numerical  scheme  based  on  adaptive  mesh  is  capable  of  handling  such  large  differences  in 
parameters  of  the  distribution  function. 


Figure  1 7.  Evolution  of  macroparameters  of  the  simulated  ion  distribution  function  in  an 
external  field  under  charge-exchange  collisions. 

3.3  Electron  isotropization  under  effect  of  elastic  collisions 

The  ratio  of  electron  to  atom  mass,  ml  MU  10”4,  is  a  small  parameter  of  the  system.  In  this 
case,  the  Lorentz  gas  model  can  be  used  to  describe  elastic  collisions  of  electrons  with  gas 
molecules.  The  collision  integral  is  a  sum  of  two  terms,  Se  =  S°e  +Sle .  The  first  term  describes 
momentum  relaxation  in  the  limit  ml M  ^  0  (Ref.  [18],  Eq.  6.1.18) 

s°e  =  4N  j  d3\S(l2  +  %  ■  l)/(£  l )  [M  +  21)  -  m\  ( 1 6) 

s2 

where  the  vector  1  corresponds  to  the  change  of  v  due  to  the  collision,  £  =  %  + 1 .  This  term 
vanishes,  S(°  =  0 ,  for  any  isotropic  distribution,  since  |  \  +  21 1=  %  .  The  differential  scattering 
cross  section  /(£,/)  is  determined  by  the  intermolecular  interaction  potentials.  For  hard 
spheres  with  diameter  a,  1(^,1)  =  a2  /  4.  The  second  term,  S'],  accounts  for  energy  exchange 
between  the  light  and  heavy  particles,  it  can  be  written  in  a  differential  (Fokker-Planck)  form. 

The  integration  (16)  takes  place  over  a  sphere  in  velocity  space  and  can  be  alternatively  written 
as 


s°e  =  JVvj  <x(£|n-ri  I )[/(£ri)-/(&n)}ra  (17) 

s2 

where  ft  is  a  velocity  angle  on  a  unit  sphere  S2  in  velocity  space,  \  =  ^ft ,  and  a  is  the 
differential  collision  cross  section,  which  depends  on  the  initial  speed  and  the  angle 
<9  =  arccos(n  •  n )  between  the  initial  and  final  electron  velocity.  When  the  scattering  is  close 
to  isotropic  cr(£,|  Cl -Cl  |)  =  cr(£) ,  the  integral  (17)  is  further  simplified 
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(18) 


seel  =  -v[f(Z,n)-lf(Z,n)dn] 

S2 

where  v  =  is  the  collision  frequency. 

This  section  describes  calculation  performed  with  the  collision  integral  (17)  for  the  simplest 
case  of  isotropic  scattering  (918).  We  first  studied  two  cases  without  external  force. 

3.3.1  Isotropization  of  an  initially  non-isotropic  distribution 

In  this  case  the  initial  distribution  function  is  has  an  elliptic  shape  with  zero  mean  velocity  and 
Tx  =  3Ty.  This  case  is  important  for  testing  how  the  implemented  collision  integral  conserves 
the  moments  of  the  distribution  function.  The  computational  grid  in  velocity  space  and  the 
initial  and  final  velocity  distribution  functions  are  shown  in  Figure  18.  The  evolution  of  the  Tx 
and  Ty  temperatures  calculated  from  the  velocity  distribution  function  is  shown  in  Figure  19. 


Figure  18.  Initial  and  final  distribution  function  for  case  with  initially  non-isotropic 

distribution. 


time*nu 

Figure  19.  Time  evolution  of  temperatures  for  the  case  with  an  initially  non-isotropic 

distribution. 

3.3.2  The  case  with  initially  non-isotropic  and  shifted  distribution  without  force 

In  this  case,  the  initial  velocity  distribution  function  has  non-zero  velocity  mean  velocity,  ^ x  = 

0.3,  and  temperature  7=0.005.  Figure  20  shows  the  dynamics  of  the  isotropization  process.  The 
corresponding  macroparameters  are  shown  in  Figure  21.  The  density  is  conserved  with  a 
machine  precision.  The  mean  velocity  drops  to  zero  in  about  two  collisional  times.  The 
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temperatures  Tx  and  Ty  evolve  towards  a  single  temperature  with  the  same  time  scale. 


Figure  20.  The  computational  grid  and  the  initial,  intermediate  and  final  distribution 
functions  (color)  under  the  action  of  the  isotropization  integral. 


time'nu  time'nu 

Figure  21.  Time  evolution  of  macroparameters  of  a  distribution  function  with  initial  non¬ 

zero  velocity. 


3.3.3  The  case  with  initially  isotropic  distribution  and  force 

Finally,  a  test  with  a  force  has  been  carried  out.  Figure  22  shows  distribution  function  at  three 
time  moments  under  the  action  of  the  isotropization  integral  (frequency  of  1)  and  in  presence 
of  force  (amplitude  of  0.25).  One  can  see  the  distribution  function  moves  in  the  direction  of  the 
force  and  at  the  same  time  evolves  in  the  c  -direction  under  the  action  of  the  isotropization 

integral  (along  |  E,  |=const  circles).  The  velocity  grid  dynamically  adapts  to  the  features  of  the 
instantaneous  distribution  function.  The  temperature  (or  average  energy)  grows  with  time  (with 
slight  difference  of  Tx  and  Ty),  as  see  in  Figure  23.  Since  there  are  no  energy  losses,  the 
temperature  will  grow  indefinitely. 
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Figure  22.  Distribution  function  at  three  time  moments  under  the  action  of  the 

isotropization  integral  and  in  presence  of  force. 


time*nu  time*nu 


Figure  23.  Time  evolution  of  macroparameters  of  a  distribution  function  under  the  action 

of  the  isotropization  integral  and  in  presence  of force. 


In  order  to  test  the  implementation  of  the  isotropization  integral,  we  have  carried  a  test  with 
force  of  a  coarse  grid.  The  results  with  Maxlevel=7  show  practically  no  difference  with  those 
on  a  fine  grid  (Maxlevel=8).  The  simulation  with  a  Maxlevel=7  took  only  30  sec  of  computing 
time. 
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4. 


FLUID  PLASMA  MODELS 


The  fluid  models  of  plasma  utilize  two  (density  and  velocity)  or  three  (density,  velocity,  and 
temperature)  equations  for  electrons  and  ions  coupled  to  electromagnetic  solvers.  During  Phase 
I,  we  advanced  fluid  plasma  models  in  two  directions. 

4.1  Basic  Plasma  Solver  with  Dynamically  Adaptive  Mesh 

We  developed  new  capabilities  for  fluid  plasma  simulations  with  dynamically  adaptive 
Cartesian  mesh  using  simple  plasma  model.  The  model  contains  transport/ionization  of 
electrons  and  ions  coupled  to  the  Poisson  equation  for  the  electric  field.  This  model  is 
appropriate  for  simulation  of  collision-dominated  high-pressure  discharges.  Our  computational 
tool  is  unique  in  several  aspects:  i)  plasma  equations  are  solved  with  dynamically  adaptive 
quadtree/octree  Cartesian  mesh,  which  provides  an  excellent  compromise  between  the 
flexibility  of  unstructured  meshes  and  the  computational  efficiency  of  structured  meshes,  ii) 
complex  boundaries  are  represented  using  the  volume-of-fluid  approach  enabling  simulations 
of  curved  electrodes,  iii)  large  dynamic  range  of  grid  refinement/coarsening  (up  to  10-12 
levels)  provides  high  resolution  of  the  streamer  fronts. 

The  efficient  and  accurate  technique  for  grid  refinement/coarsening  is  very  important  for 
simulations  of  streamers  with  dynamically  adaptive  mesh.  It  was  found  empirically  that  for 
reliable  results,  the  following  criteria  had  to  be  used: 

\neM)-neM^\  <  0*05max(«eJ)  and  nej(rj)/nej(rj+l)  <  2 

|E(ry)|/|E(ry+1)|  <  2  and  v^/v^)  <  5 

where  ne  and  nt  are  the  electron  and  ion  densities,  E  is  the  electric  field  and  v.  is  the  ionization 
frequency,  and r.  and  rj+l  denote  neighboring  cells. 

An  example  of  3D  simulations  of  streamer  development  from  an  initial  perturbation  is  shown 
in  Figure  24a.  No  branching  of  streamers  has  been  observed  in  our  simulations.  No  streamers 
develop  in  high  electric  fields,  dumb-bell  shape  structures  appear  instead.  Figure  24b  shows 
2D  axi-symmetric  streamer  development  near  a  high-voltage  needle-like  cathode.  The  streamer 
propagates  with  a  velocity  by  an  order  of  magnitude  higher  than  the  electron  drift  velocity  in 
the  highest  electric  field  in  the  gap.  The  calculated  velocity  agrees  with  classical  estimates  for 
the  fast  streamer  velocity.  Figure  24c  illustrates  streamer  development  from  two  spatially 
separated  avalanches.  Both  negative  and  positive  streamers  expand  with  velocities  higher  than 
the  electron  drift  velocity.  The  negative  streamer  velocity  is  3  times  higher  than  the  drift 
velocity,  the  positive  streamer  velocity  is  about  2  drift  velocities.  The  neighbor  streamers  get 
closer  during  their  development;  the  field  increases  in  the  gap  between  the  streamer  tips,  and 
ionization  takes  place  along  the  line  connecting  the  positive  and  negative  tips. 

In  all  cases,  the  calculated  fields  near  the  streamer  tips  were  of  the  order  of  the  critical  field 
corresponding  to  saturation  of  the  ionization  frequency  with  respect  to  the  field  strength.  This 
implies  that  non-local  ionization  should  play  an  important  role  in  these  regions. 
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Figure  24.  a)  3D  simulation  of  a  streamer  development  from  an  initial  perturbation  in  N2: 

p=400  Torr,  E=34  kV/cm,  b)  2D  axi-symmetric  simulation  of  streamer  development  between  a 
needle-like  elliptic  cathode  (size  2x1  cm)  and  a  flat  anode:  p=7 60  Torr,  cathode  potential  -600 
kV,  c)  interaction  of  two  avalanches  during  their  transition  to  streamers:  p=7 60  Torr,  E=65 
kV/cm,  computational  domain  2x1  cm  and  initial  plasma  density  5xl010  cm3. 


4.2  Adding  Forces  in  the  gas  kinetic  continuum  solvers 

We  have  studied  numerical  algorithms  associated  with  external  forces  in  the  gas  kinetic  Euler 
and  Navier-Stokes  solvers.  The  implementation  of  forces  in  gas  kinetic  schemes  and 
difficulties  associated  with  “numerical  heating”  are  described  in  Ref.  [19].  For  tackling  this 
effect,  we  used  the  approach  proposed  in  Ref.  [20]. 

In  order  to  test  different  algorithms,  we  have  considered  a  problem  of  a  stationary  gas  in 

T  17TX 

periodic  field  with  potential  (p  =  -% — sin( - ).  The  analytical  solution  to  this  problem 

In  L 

corresponds  to  zero  mean  velocity,  constant  gas  temperature,  and  the  density  distributed 
according  to  the  Maxwell-Boltzmann  relation  {p  ~  exp(-^/T)).  The  most  difficult  part  in  this 
benchmark  test  is  to  achieve  small  velocity  and  constant  temperature.  We  have  used  3  different 
algorithms  to  account  for  the  force  terms  in  our  NS  solver:  the  first  approach  uses  a  limiter 
(Van-Leer),  the  second  case  is  without  a  limiter  and  the  third  one  is  from  Ref.  [19]),  which 
accounts  for  the  terms  with  the  force  inside  the  distribution  function.  The  results  using  these 


22 


8840/Final 


methods  are  compared  in  Figure  25.  One  can  see  that  for  the  third  method,  the  smallest 
velocity  values  could  be  obtained  and  the  temperature  is  close  to  a  constant.  The  results  are 
very  encouraging  since  they  show  that  it  is  possible  to  achieve  high  precision  by  using  the 
appropriate  numerical  implementation  for  the  force  term  in  gas  kinetic  schemes.  This  approach 
is  planned  for  future  coupling  of  the  Boltzmann  solution  with  the  NS  solution  in  the  presence 
of  forces  in  realistic  plasma  systems. 
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5. 


CONCLUSION 


During  Phase  I,  we  have  investigated  methods  of  Renormalization  Group  and  Slow  Invariant 
Manifolds  in  application  to  plasma  modeling.  Using  the  symmetry  group  for  the  Vlasov- 
Maxwell  system,  we  obtained  exact  solutions  for  the  dynamics  of  boundary  layer  in 
collisionless  plasma.  The  interpretation  of  the  solutions  in  terms  of  invariants  of  the  group 
transformation  and  invariant  manifolds  was  achieved  -  the  particle  distribution  functions  were 
related  to  invariants  of  the  group  transformation.  The  road  towards  obtaining  approximate 
symmetries  due  to  small  parameter  (electron/ion  mass  ratio)  in  the  plasma  was  identified.  We 
prepared  detailed  work  plan  for  further  development  and  application  of  this  methodology  for 
weakly  ionized  collisional  plasmas  in  Phase  II.  We  have  analyzed  and  tested  several  numerical 
methods  of  solving  the  Vlasov  and  Boltzmann  kinetic  equations  for  charged  particles  in 
weakly  ionized  plasma.  We  found  that  cubic  spline  interpolation  provides  superior  results 
compared  to  the  second  order  scheme  with  van  Leer  limiter.  The  Eulerian  Vlasov  solver  with 
cubit  interpolation  reproduced  analytical  results  for  the  shape  of  the  velocity  distribution 
function  with  high  accuracy  using  a  reasonable  number  of  mesh  points.  We  have  evaluated 
numerical  methods  of  solving  hydrodynamic  plasma  equations  for  different  types  of 
computational  grid.  We  have  developed  a  minimal  plasma  solver  with  dynamically  adaptive 
mesh  capabilities  and  demonstrated  this  solver  for  the  problem  of  streamer  development  in 
different  electric  fields.  We  have  submitted  a  paper  to  IEEE  Transactions  on  Plasma  Science 
entitled  "Streamer  Simulations  with  Dynamically  Adaptive  Cartesian  Mesh",  which  was 
accepted  for  publication. 
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